# Read in Data
all_mfri = list.files('./Wildfire_MFRI/',pattern = '.tif',full.names = T)
for(file in all_mfri){
object_name = file_path_sans_ext(basename(file))
assign(object_name, raster(file))
}
CC4a_reg = read_sf('./Boundries/CC4a_RegionsSub.shp')
CC4a_reg = st_transform(CC4a_reg, "+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
# dissolve multipart Sierra Mountain feature
CC4a_reg = ms_dissolve(CC4a_reg,field = 'Region')
# Get stats for 01-25 & 26-50 MFRIs cap at 500 years
all_2001_2025 = stack(all_mfri[grepl('2001_2025',all_mfri)])
all_2026_2050 = stack(all_mfri[grepl('2026_2050',all_mfri)])
# years for maximum plotted MFRI
capat=500
# calculate stats for stacks and cap MFRI
# writes out files with 3 different postfix _stat for statisic (mean etc), _df for dataframe for ggplot, _capped for raster capted at value capat
summary_functions = c('min','max','mean' )
for(summary in summary_functions){
for(stac in c('all_2001_2025','all_2026_2050')){
assign(paste(summary,stac,sep='_'),do.call(summary,list(x=get(stac),na.rm=T)))
capper = get(paste(summary,stac,sep='_'))
# write out df for ggplot
capper.df = data.frame(rasterToPoints(capper))
names(capper.df) =c("lon", "lat","MFRI")
assign(paste(summary,stac,'df',sep='_') ,capper.df)
# cap at catat yrs
capper[capper>capat]=capat
assign(paste(summary,stac,'capped',sep='_'),capper)
}}
# write out df for ggplot
MFRI_76_00.df = data.frame(rasterToPoints(MFRI_76_00))
names(MFRI_76_00.df) =c("lon", "lat","MFRI")
all_2001_2025[all_2001_2025>500]=500
all_2026_2050[all_2026_2050>500]=500
plot(all_2001_2025)
plot(all_2026_2050)
mean_chg_76_25 = MFRI_76_00.df
mean_chg_76_25$MFRI = mean_all_2001_2025_df$MFRI-mean_chg_76_25$MFRI
mean_chg_76_25$MFRI[ mean_chg_76_25$MFRI >350]=350
rng= range(mean_chg_76_25$MFRI)
ggplot()+geom_raster(data=mean_chg_76_25,aes(x=lon,y=lat,fill=MFRI))+
scale_fill_gradientn(colours= c("#cc0000", "#cc0000" , 'grey', "#339933","#339933" ), #colors in
limits=c(-350, 350))+ geom_sf(data=CC4a_reg,colour = "grey30", fill = NA,size=.75) + #same limits for plots
ggtitle('Change in MFRIs 2000 - 2025 mean model run')+coord_sf()
for(region_aoi in CC4a_reg$Region){
aoi = st_bbox(CC4a_reg[CC4a_reg$Region==region_aoi,])
aplot = ggplot()+ geom_raster(data=mean_chg_76_25,aes(x=lon,y=lat,fill=MFRI))+
scale_fill_gradientn(colours= c("#cc0000", "#cc0000" , 'grey', "#339933","#339933" ),
limits=c(-350, 350))+ geom_sf(data=CC4a_reg,colour = "grey30", fill = NA,size=.75) +
ggtitle(paste(region_aoi,'\nChange in MFRIs 2000 - 2025 \nMean model run'))+
coord_sf( xlim=c(aoi[1],aoi[3]),ylim=c(aoi[2],aoi[4]) )
plot(aplot)
}
mean_chg_76_50 = MFRI_76_00.df
mean_chg_76_50$MFRI = mean_all_2026_2050_df$MFRI-mean_chg_76_50$MFRI
mean_chg_76_50$MFRI[ mean_chg_76_50$MFRI >350]=350
rng= range(mean_chg_76_50$MFRI)
ggplot()+geom_raster(data=mean_chg_76_50,aes(x=lon,y=lat,fill=MFRI))+
scale_fill_gradientn(colours= c("#cc0000", "#cc0000" , 'grey', "#339933","#339933" ),
limits=c(-350, 350))+ geom_sf(data=CC4a_reg,colour = "grey30", fill = NA,size=.75)+
ggtitle('Change in MFRIs 2000 - 2050 mean model run')+coord_sf()
for(region_aoi in CC4a_reg$Region){
aoi = st_bbox(CC4a_reg[CC4a_reg$Region==region_aoi,])
aplot = ggplot()+ geom_raster(data=mean_chg_76_50,aes(x=lon,y=lat,fill=MFRI))+
scale_fill_gradientn(colours= c("#cc0000", "#cc0000" , 'grey', "#339933","#339933" ),
limits=c(-350, 350))+ geom_sf(data=CC4a_reg,colour = "grey30", fill = NA,size=.75) +
ggtitle(paste(region_aoi,'\nChange in MFRIs 2000 - 2050 \nMean model run'))+
coord_sf( xlim=c(aoi[1],aoi[3]),ylim=c(aoi[2],aoi[4]) )
plot(aplot)
}
the_fun = median # function used to summarize raster values by polygons
region_code = data.frame(ID=seq(1,9),as(CC4a_reg,'Spatial')@data$Region)
full_mean_stack = stack(MFRI_76_00,mean_all_2001_2025,mean_all_2026_2050)
names(full_mean_stack) =c('2000','2025','2050')
extract_full_mean_df = extract(full_mean_stack, as(CC4a_reg,'Spatial'), fun=the_fun, na.rm=T, df=T)
extract_full_mean_df = left_join(region_code,extract_full_mean_df,by='ID') %>% select(-ID)%>%melt()
extract_full_mean_df$Year = as.numeric(substr(as.character(extract_full_mean_df$variable),2,5))
names(extract_full_mean_df)=c('Region','variable','value','Year')
the_fun = min # function used to summarize raster values by polygons
full_min_stack = stack(MFRI_76_00,min_all_2001_2025,min_all_2026_2050)
names(full_min_stack) =c('2000','2025','2050')
extract_full_min_df = extract(full_min_stack, as(CC4a_reg,'Spatial'), fun=the_fun, na.rm=T, df=T)
extract_full_min_df = left_join(region_code,extract_full_min_df,by='ID') %>% select(-ID)%>%melt()
extract_full_min_df$Year = as.numeric(substr(as.character(extract_full_min_df$variable),2,5))
names(extract_full_min_df)=c('Region','variable','value','Year')
ggplot()+geom_smooth(data=subset(extract_full_min_df, Region != 'Inland South' ),aes(x=Year,y=value,color=Region,group=Region),size=1.25,alpha=.6) +ggtitle('Minimum observed MFRI of Min(all_models) by region - \n omitting Inland South')
ggplot()+geom_smooth(data=subset(extract_full_mean_df, Region != 'Inland South' ),aes(x=Year,y=value,color=Region,group=Region),size=1.25,alpha=.6) +ggtitle('Median observed MFRI of Mean(all_models) by region - \n omitting Inland South')